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Abstract 

The Isthmus of Tehuantepec has played an important role in shaping the avian 
diversity of Mexico, as well as the rest of the Western Hemisphere. It has been 
both a barrier and a land connector between North and South America for 
many groups of birds. Furthermore, climatic change over the Pleistocene has 
resulted in ecological fluctuations that led to periods of connection and isola- 
tion of the highlands in this area. Here we studied the divergence of two species 
of orioles whose distribution in the highlands is separated by the lowlands of 
the Isthmus of Tehuantepec: Icterus graduacauda (west of the Isthmus) and 
Icterus chrysater (east of the Isthmus). We sequenced multiple loci (one mito- 
chondrial gene and six nuclear introns) and performed coalescent analyses 
(Isolation with Migration) to test whether their divergence resulted from prior 
occupancy of the ancestral area followed by a vicariant event or recent dispersal 
from one side or the other of this Isthmus. Results strongly indicate a vicariant 
event roughly 300,000 years ago in the Pleistocene followed by little or no gene 
flow. Both mitochondrial and nuclear genes show that the Isthmus of Tehuan- 
tepec is a strong barrier to gene flow. Thus, these two species appear to not 
exchange genes despite their recent divergence and the close geographic prox- 
imity of their ranges. 



doi: 10.1002/ece3.768 



Introduction 

The Great American Biotic Interchange, which began 
around 3.5 million years ago with the completion of the 
Panama land bridge, resulted in numerous plants and 
animals moving between North and South America (Weir 
et al. 2009). Thus, these previously isolated fauna from 
temperate and tropical regions were intermixed, resulting 
in unique evolutionary and ecological consequences from 
these regions. Even though the Isthmus of Panama is the 
best-known portion of the land bridge, another lowland 
area, the Isthmus of Tehuantepec in southern Mexico, is 
also an important part of the connection between the 
Americas, the Isthmus of Tehuantepec in southern 
Mexico. 



The Isthmus of Tehuantepec (Fig. 2) is a strip of land 
that constricts to a width of roughly 200 km between 
the Atlantic and the Pacific. Furthermore, the Oaxacan 
highlands to the west and the highlands of Chiapas to the 
east (Binford 1989; Barrier et al. 1998) are only separated 
by 50 km. Elevation drops from 2000+ m in the high- 
lands on either side of the isthmus to 200 m in the low- 
lands, making it an important barrier to dispersal for 
species of higher elevation (Sullivan et al. 2000; Castoe 
et al. 2009; among others). The Isthmus of Tehuantepec 
provides a model location for studying modes of specia- 
tion as it can serve as a barrier for some species, while 
acting as a connector for others. 

The two fundamental modes of geographic speciation 
are vicariance and dispersal (founder event speciation). 
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The Isthmus of Tehuantepec may have fragmented two 
populations from a widespread ancestral population and 
isolated them for enough time that they diverge into dis- 
tinctive lineages as a result of the vicariant event. On the 
other hand, speciation by dispersal could occur with the 
isthmus acting as a connector, if some individuals from 
either side moved across the barrier and founded a new 
geographically isolated population. 

Following either speciation mode, the allopatric popu- 
lations could experience little or no gene flow and 
become increasingly reproductively isolated (Coyne and 
Orr 2004). However, a growing body of research demon- 
strates that even after a population splits, gene flow can 
still occur in time periods following initial divergence 
(Hey 2005). For example, climatic oscillations (e.g., in 
temperature, precipitation, and sea levels) during the 
Pliocene or Pleistocene likely caused the Isthmus of Tehu- 
antepec to experience alternating periods of connection 
and disconnection between different highland habitats 
(Duellman 1960; Zarza et al. 2008). 

A number of studies have examined the role the Isth- 
mus of Tehuantepec has played in the speciation of ter- 
restrial species, including birds. The Isthmus of 
Tehuantepec likely acted as a barrier to dispersal, separat- 
ing populations to the west and east, for montane birds 
such as Chlorospingus ophthalmicus (Sanchez-Gonzalez 
et al. 2007), Lampornis amethystinus (Cortes-Rodriguez 
et al. 2008), and Picoides villosus (Klicka et al. 2011). 
Another study based on mitochondrial DNA examined 
several bird species with distributions on either side of 
the Isthmus and concluded that fluctuations in the cli- 
mate of montane habitats had fractured the bird fauna 
during the late Pliocene and Pleistocene (Barber and 
Klicka 2010). However, other studies have found that the 
Isthmus of Tehuantepec has not acted as a barrier for dis- 
persal given the lack of genetic differentiation between 
populations of other highland bird species such as Lepi- 
docolaptes affinis (Arbelaez-Cortes et al. 2010) and Buarre- 
mon species (Navarro-Sigiienza et al. 2008) or lowland 
species such as Icterus pustulatus. 



Traditionally, this type of biogeographic research, 
including both species-level phylogenetics and within- 
species phylogeography, has used mitochondrial DNA as 
a way to understand geographic variation within and 
between closely related species. However, mtDNA alone is 
often insufficient to elucidate phylogeographic patterns 
(Funk and Omland 2003; Peters et al. 2012a). The evolu- 
tionary history of a gene might not be congruent with the 
evolutionary history of a species due to lineage sorting of 
ancestral polymorphisms or introgressive hybridization; 
thus, relying on mitochondrial DNA alone can lead to 
incorrect conclusions about the evolutionary history of 
lineages (Moore 1995; Jacobsen and Omland 2012). Use 
of multiple independent loci provides a more complete 
picture of population histories because mutation, genetic 
drift, and natural selection operate independently in 
unlinked loci (Peterson and Nyari 2008; Peters et al. 
2012b). 

One current challenge is to study the relationship 
between shared genetic variation, the timing of diver- 
gence, and the geographic context of divergence along 
with other processes such as ongoing gene flow (Won 
and Hey 2005). Coalescent methods can help answer 
questions at the level of population divergence; indeed, 
these methods are best applied to within species or 
between recently diverged taxa (Rosenberg and Nordborg 
2002). The program Isolation with Migration (IM; Hey 
and Nielsen 2004) uses data from unlinked genes to 
obtain posterior probability distributions for several 
demographic parameters including effective population 
size, migration, and time since divergence, and together, 
they capture many of the dynamics of the early stages of 
population differentiation (Hey and Machado 2003). 

Figure 1 illustrates Icterus graduacauda (A) and Icterus 
chrysater (B), two sister species of tropical orioles whose 
distributions are separated by the Isthmus of Tehuante- 
pec, and with very different colorations. The two species 
are roughly 1.6% divergent in mitochondrial cytochrome 
b and ND2 (Omland et al. 1999), suggesting their split 
occurred less than one million years ago (Fleischer et al. 




(photographed 
js chrysater 



4378 



© 201 3 The Authors. Ecology and Evolution published by John Wiley & Sons ttd. 



N. Cortes-Rodriguez ef al. 



No migration between two closely related Orioles 



1998). These two orioles thus seem to be more closely 
related than any other sister species in the genus (Omland 
et al. 1999; Jacobsen et al. 2010). As other orioles more 
deeply divergent hybridize extensively, at least some gene 
flow between these parapatric taxa seems likely. Further- 
more, these two orioles are among the few species within 
the genus Icterus that inhabit highland habitats, reaching 
as high as 3000 m (Jaramillo and Burke 1999). Icterus 
graduacauda inhabit highlands west of the Isthmus of 
Tehuantepec from Texas to Oaxaca. Icterus chrysater 
inhabits highlands east of the Isthmus of Tehuantepec 
from eastern Chiapas to northern South America (Howell 
and Webb 1995). Even though these two orioles inhabit a 
wide variety of vegetation types along their range, near the 
edges of their distribution (close to the Isthmus of Tehu- 
antepec) I. graduacauda inhabits humid pine-oak forests 
and cloud forests above 2500 m, whereas Icterus chrysater 
inhabits arid and semiarid pine-oak forests above 1500 m 
(Binford 1989; Jaramillo and Burke 1999). There are no 
known sympatric populations of the two species, and no 
hybridization has been documented or suggested. 

In this study, we used IM to test two alternative 
hypotheses. HI) the two species resulted from vicariant 
speciation where the Isthmus of Tehuantepec as a barrier 
split the widespread ancestral population into two daugh- 
ter populations. H2) the two species resulted from a 
founder event, in which the ancestral population was 
located on one side of the Isthmus and a small number of 
individuals crossed the barrier, founding a new daughter 
population. The direction of colonization could either 
have been from west to east or from east to west. Follow- 
ing Hey (2005) and Peters et al. (2008), we used the split- 
ting parameter s to estimate the percentage of the 



ancestral population that contributed to each daughter 
population. 

Methods 

Taxon sampling and laboratory procedures 

We obtained frozen tissues and toepads from individuals 
of Icterus graduacauda (N = 38) and Icterus chrysater 
{N = 49) collected from across each species' range 
(Fig. 2). DNA was extracted from both fresh muscle and 
toepad samples using the DNeasy Blood and Tissue Kit 
(Qiagen, Valencia, CA). To improve DNA yield from toe- 
pad samples (age of the samples ranged between 1940 
and 1980), we added 20 fiL of dithiothreitol (DTT) to the 
sample in addition to 20 /jL of proteinase K. 

We sequenced a total of seven loci, including 344 bp of 
the mitochondrial control region (CR-domain I) for all 
individuals and six intronic regions of nuclear DNA for a 
subset of individuals (a-ENO, RDP2, GAPDH, TGFB2, 
MUSK, and SLC9, see Table 1 for details). PCR was car- 
ried out on the Gene Amp PCR System 9700 (Applied 
Biosystems, Foster City, CA), and products were verified 
on a 1% agarose gel, visualized by staining with ethidium 
bromide, and cleaned using QIAquick PCR purification 
kits (Qiagen). We sequenced these products using ABI's 
BigDye v.3 Terminator Cycle Sequencing Kit on the Gene 
Amp PCR system 2400 (Applied Biosystems, Foster 
City, CA). Excess dye terminators were removed by the 
EDTA-ethanol precipitation protocol as recommended by 
Applied Biosystems, and sequences were run out at 
University of Maryland, Baltimore County (UMBC) on 
an ABI 3100 DNA sequencer (Applied Biosystems). All 
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Table 1. Characteristics of loci used in this study. One mitochondrial locus and six nuclear introns with the total number of base pairs, chromo- 
some location, annealing temperature, and the primer source. 



Locus 


Total bp 


Chromosome 


Annealing temperature 


Primer source 


Control region (CR) 


344 


Mitochondria 


52°C 


Kondo et al. (2004) 


GADPH 


283 


Chromosome 1 


55°C 


Primmer et al. (2002) 


TGFB 


544 


Chromosome 3 


65°C 


Bures et al. (2002) 


RDP2 


297 


Chromosome 12 


56°C 


Waltari and Edwards (2002) 


ENO 


252 


Chromosome 21 


60°C 


Kondo et al. (2008) 


MUSK 


467 


z-chromosome 


50°C 


Clark and Witt (2006) 


SLC9 


434 


z-chromosome 


56°C 


Backstrom et al. (2006) 



sequences were edited and aligned using Sequencher 4.1 
(Genecodes Corporation, Ann Arbor, MI). 

Population structure and selection 

We constructed haplotype networks for all loci using the 
median-joining algorithm in Network 4.5.1.6 (Bandelt 
et al. 1999). For each species, we calculated the following: 
nucleotide diversity (n), nucleotide polymorphisms (9), 
Fu's Fs (Fu 1997), R 2 (Ramos-Onsins and Rozas 2002), 
Tajima's D (Tajima 1989), and mismatch distributions 
(with their associated Tau and Theta) using the program 
DnaSP 4.0 (Rozas et al. 2003). All of the above calcula- 
tions were tested; however, Tajima's D and R 2 are the rec- 
ommended statistics (Ramirez-Soriano et al. 2008). We 
then estimated the level of differentiation between species 
and populations within each species for each locus (for 
within- and between-species variation) using analysis of 
molecular variance (AMOVA) in Arlequin v3.01 (Excof- 
fier et al. 2006). Populations within each species were 
established based on sample localities. We used calcula- 
tions in AMOVA that took into account the number of 
mutations between haplotypes (^-statistics). In order to 
calculate the allelic richness for each nuclear locus, we 
standardized the sample sizes to the smallest sample using 
the software Rarefaction calculator (University of Alberta, 
Alberta, Canada). 

Coalescent analyses 

We used coalescent analysis of six unlinked nuclear loci 
to study the divergence history of the sister species 
I. graduacauda and I. chrysater (the control region of the 
mtDNA was not included). The program Isolation with 
Migration (IM; Hey and Nielsen 2004) jointly estimates 
demographic parameters to obtain posterior probability 
distributions from unlinked genes (Hey 2007). These 
parameters offer the opportunity of capturing the dynam- 
ics of a population during the early stages of differentia- 
tion (Hey 2005). The model estimates the following 
demographic parameters scaled to the per-locus mutation 



rate: effective population sizes (6i = population 1, 
6 2 = population 2, and 9 A = ancestral population), migra- 
tion rates (m t and m 2 ), and the time of population split- 
ting (r). IM can also estimate the fraction of the ancestral 
population that founded each population (s and 1-s). The 
program IM operates under the following assumptions: 
(1) loci are neutral, (2) no recombination within loci, (3) 
free recombination between loci, (4) panmictic ancestral 
and descendant populations, and (5) constant migration 
rates over time (Hey and Nielsen 2004). 

Allelic states at nuclear loci were inferred using PHASE 
v.2.1 (Stephens et al. 2001; Stephens and Donnelly 2003). 
After being phased, every locus was tested for evidence of 
recombination using IMgc (Woerner et al. 2007). IMgc 
infers recombination from violations of the four gamete 
rule within a sequence alignment by testing compatibility 
between pairs of polymorphic sites (Hudson and Kaplan 
1985). Inheritance scalars were defined depending on 
the mode of inheritance for each locus (autosomal = 1, 
and Z-linked = 0.75), and the infinite-sites mutation 
model was used. To improve mixing of the Markov 
chains (to facilitate convergence), we ran multiple heated 
chains and kept monitoring the autocorrelation and esti- 
mates of ESS (effective sample sizes; Hey 2007). We used 
Isolation with Migration Analytical (IMa; Hey and Niel- 
sen 2007) to obtain migration, effective population size 
(ancestral population and daughter populations) and 
divergence time. Because IMa does not calculate the split- 
ting parameter s, we also used Isolation with Migration 
(IM; Hey and Nielsen 2004). 

We repeated the analysis twice using the same priors 
but different seeds in each one of the runs. The runs 
included 1,000,000 generations, using 30 heated chains 
and sampling every 100 generations. To convert the 
estimated parameter into biologically informative values, 
we need the generation time and mutation rates. We 
assumed an average generation time for passerine birds of 
1.7 years (Sasther et al. 2005). We used a mutation rate 
average of 3.6 x 10~ 9 s/s/y (Axelsson et al. 2004). These 
mutational rates were provided in the IM input file and 
used to convert the parameter estimates into biological 
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Figure 3. Haplotype networks of the six nuclear loci sequenced. Note that all loci show haplotype sharing between species. 



informative values. We ran the model for 100,000,000 
generations. 

Results 

Population structure, selection, and 
recombination 

Haplotype networks show allele sharing between 7. gradu- 
acauda and I. chrysater across all nuclear loci (Fig. 3). 
Nearly all loci display a similar distinctive star-shaped 
pattern, with one or two shared central ancestral haplo- 
types surrounded by several rare haplotypes that are not 
shared between species, which could be the result of a 
recent population expansion or positive selection (Schnei- 
der and Excoffier 1999). Also apparent in the haplotype 
networks, all introns and the control region show high 
levels of haplotype diversity (except RDP2, which showed 
only a few haplotypes per species). Haplotype diversity, 
nucleotide diversity, and allelic richness vary between loci 
in both I. graduacauda and I. chrysater (Table 2). 

The control region haplotype network (Fig. 4) shows 
no haplotype sharing between I, graduacauda and 
I. chrysater, which could be due to the rapid sorting rate 
of mtDNA compared with nDNA. These two groups are 



on average 2.6% divergent. Assuming a rough divergence 
rate of 4.0% per million years for control region (Flei- 
scher et al. 1998), the mtDNA of these two species started 
diverging around 600,000 years ago. I. graduacauda had 
substantially higher mitochondrial haplotype and nucleo- 
tide diversity than I. chrysater (haplotype diversity: 0.86 
vs. 0.46; nucleotide diversity: 0.0060 vs. 0.0015; see 
Table 2). However, nucleotide diversity between the two 
species varies greatly between the introns; for example, 
I. graduacauda has a higher nuclear diversity in GADPH 
and SLC9, whereas I. chrysater has a higher nuclear diver- 
sity in a-ENO, RDP2, TGFB2, and MUSK (Table 2). 

The AMOVA results (Table 3) indicate high genetic 
variation within populations of the two species (I. gradu- 
acauda and I. chrysater) for most of the nuclear introns 
and with low $ st values (except for GADPH). As 
expected, given fixed mtDNA differences, genetic varia- 
tion between species is more apparent when mtDNA is 
analyzed, with species explaining 84.9% of the total varia- 
tion and a $ st of 0.88, which is expected given the lack of 
haplotype sharing between these two oriole species, as can 
be observed in the haplotype network. Although sample 
sizes are generally too low to warrant AMOVA on 
individual loci, all nuclear loci show low amounts of vari- 
ation explained by the between-species comparison. 
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Table 2. Estimates of haplotype diversity (including number of haplotypes), nucleotide diversity, and allelic richness (standardized to the smallest 
sample size for each locus) for the seven loci sequenced for Icterus graduacauda and Icterus chrysater. The highest value for nucleotide diversity 
each locus is shown in bold. 



Locus 


No. of 

individuals 

(/. graduacauda) 


No. of 
individuals 
(/. chrysater) 


Haplotype 
diversity 

(/. graduacauda) 


Haplotype 
diversity 
(/. chrysater) 


Nucleotide 
diversity 

(/. graduacauda) 


Nucleotide 
diversity 
(/. chrysater) 


Allelic richness 
(/. graduacauda) 


Allelic richness 
(/. chrysater) 


CR 


30 


26 


0.86 (15 hap) 


0.46 (4 hap) 


0.0060 


0.0015 


13.38 ± 0.92 


4 ± 0.00 


RDP2 


26 


24 


0.14 (2 hap) 


0.61 (6 hap) 


0.0005 


0.0026 


1.99 ± 0.05 


6 ± 0.00 


SLC9 


24 


18 


0.77 (8 hap) 


0.57 (7 hap) 


0.0029 


0.0022 


10.04 ± 1.32 


7 ± 0.00 


MUSK 


27 


23 


0.72 (7 hap) 


0.93 (12 hap) 


0.0026 


0.0048 


6.63 ± 0.54 


12 ± 0.00 


TGFB 


28 


26 


0.65 (6 hap) 


0.91 (12 hap) 


0.0015 


0.0038 


4.85 ± 0.36 


12 ± 0.00 


ENO 


60 


50 


0.74 (8 hap) 


0.85 (16 hap) 


0.0044 


0.0069 


7.33 ± 0.73 


16 ± 0.00 


GADPH 


42 


24 


0.73 (11 hap) 


0.69 (7 hap) 


0.0036 


0.0039 


9.43 ± 1.28 


7 ± 0.00 




Figure 4. Haplotype network for the mitochondrial control region. 
Circles are drawn proportional to number of haplotypes (note: no 
haplotype sharing is observed). 

Table 4 summarizes the results of the three tests of 
neutrality. Most loci show negative Tajima's D values 
(with the exception of a-ENO for /. graduacauda and 
GAPDH for I. chrysater). Also, Fu's Fs values are negative 
in both species for all loci. However, none of the loci 
used in this study depart significantly from neutrality 
(P > 0.05). 



Coalescent analysis using Isolation with 
Migration 

We found some evidence of recombination at all loci 
(except mtDNA) that were filtered using IMgc online. IM 
(Hey and Nielsen 2004) analysis of the six nuclear intron 
loci strongly indicates little to no gene flow between I. grad- 
uacauda and I. chrysater, with a peak estimate of migration 
at zero migrants per generation in both directions across 
the Isthmus of Tehuantepec (Fig. 5; I. graduacauda 95% 
HPD 0-0.002 and I. chrysater 95% HPD 0-0.004). The 
effective population size of I. graduacauda is about 300,000 
individuals (95% HPD 150,000-620,000), whereas Ne for 
I. chrysater is about 600,000 individuals (95% HPD 
320,000-1,200,000). The ancestral population size was esti- 
mated as 180,000 individuals (95% HPD 45,700-256,000) 
individuals. Also, the multilocus IM analysis suggests that 
I. graduacauda and I. chrysater split roughly 300,000 years 
ago (95%HPD 176,000-860,000). This time estimate indi- 
cates that these two species likely diverged sometime during 
the Pleistocene when climate change and vegetation shifts 
due to glacial cycles affected many other species as well 
(Hewitt 2004). Moreover, the splitting parameter s shows a 
peak at 0.84 signifying that 84% of the ancestral population 
contributed to I. chrysater, while 16% of the ancestral pop- 
ulation contributed to I. graduacauda (Fig. 6). 



Table 3. AMOVA results and <J> st value significance and percentage of variation explained between species. 









% of variation 


% of variation among 


% of variation 


Locus 


4>st 


P-value 


between species 


populations within groups 


within populations 


CR 


0.88 


0.00000 ± 0.00000 


84.91 


3.38 


11.71 


RDP2 


-0.068 


0.96090 ± 0.00529 


13.47 


-20.31 


106.84 


SLC9 


0.36 


0.00293 ± 0.00216 


4.87 


31.51 


63.63 


MUSK 


0.42 


0.00000 ± 0.00000 


3.69 


38.48 


57.83 


TGFB 


0.23 


0.00000 ± 0.00000 


2.36 


21.08 


76.56 


ENO 


0.14 


0.00391 ± 0.00185 


0.83 


13.38 


85.79 


GADPH 


0.60 


0.00000 ± 0.00000 


-6.00 


66.35 


39.53 
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Table 4. Test for deviation from selective neutrality and constant population sizes between Icterus graduacauda and Icterus chrysater. 



Locus 


Tajima's D 




Fu's Fs 




R2 




/. graduacauda 


/. chrysater 


/. graduacauda 


/. chrysater 


/. graduacauda 


/. chrysater 


CR 


-0.249 


-1.04 


-9.727* 


-1.534** 


0.115 


0.092 


RDP2 


-0.713 


-1.243 


-0.317 


-2.813* 


0.073 


0.086 


SLC9 


-0.739 


-1.791 


-3.66** 


-3.955* 


0.098 


0.085 


MUSK 


-0.248 


-0.017 


-1.705 


-5.003* 


0.122 


0.129 


TGFB 


-0.751 


-0.506 


-1.524 


-5.538 


0.06 


0.106 


ENO 


0.096 


-0.363 


-2.285 


_g 75** 


0.11 


0.093 


GADPH 


-1.006 


0.176 


-6.175** 


-1.928 


0.073 


0.139 



Significance of the observed values is indicated at 0.05 > 0.01 P (*) and P > 0.001 (**). Tajima's D and Fu's F show that the sequences are 
negative for the seven loci in both species. 



Figure 5. Migration probabilities from Icterus 
chrysater into Icterus graduacauda (gray line). 
And migration Icterus graduacauda into Icterus 
chrysater (black line) and migration (note: 
highest posterior probability for both directions 
is zero migrants per generation). 
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Figure 6. Posterior distribution of the s 
splitting parameter from the Isolation with 
Migration model based on the six nuDNA 
introns. 



Discussion 

Vicariance across the Isthmus of 
Tehuantepec and no gene flow between 
Icterus graduacauda and Icterus chrysater 

Multilocus coalescent methods can be used to test 
hypotheses regarding the signatures of vicariance versus 
founder events across geographical barriers. Also, IM 



0.8 0.6 0.4 0.2 0 

(splitting parameter) % of ancestral population 
that contributed to Icterus graduacauda 

results suggest high Ne estimates for both species 
(I. graduacauda = 300,000 and I. chrysater = 600,000), 
and an s parameter of 16% suggests that a reasonably 
large fraction of the ancestral population contributed to 
both species at the time of the split. These results thus 
suggest that a widespread ancestral species was divided by 
the Isthmus of Tehuantepec into two daughter popula- 
tions. Although there is no absolute threshold value for a 
founder event, when using the parameter s, 5% is a 
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reasonable rule of thumb; furthermore, our finding of 
16% is considerably higher than the values observed in 
other studies in which researchers found evidence for a 
founder event scenario. For example, Hey (2005) showed 
that the new world population of humans was found 
only by less than 1% (s = 0.99) from the Asian popula- 
tion. Similarly, Peters et al. (2008) found that the North 
American populations of Gadwalls (Anas strepera) were 
colonized by less than 1% (s = 0.04) of the Eurasian 
population. 

Surprisingly, despite the narrowness of the strip of 
lowland separating the distributions of I. graduacauda 
and I. chrysater, we found no evidence of interspecific 
gene flow across the Isthmus of Tehuantepec. Gene flow 
between the two species would have seemed likely 
because the distance between our samples from the Sierra 
Madre Oriental (I graduacauda) and the Sierra Madre 
del Sur (I. chrysater) is only about 50 km. However, dis- 
persal across this valley seems rare or nonexistent based 
on IM estimates of <0.005 migrants per generation from 
either side; if any individuals do cross, they may be 
unable to establish territories and/or obtain mates. These 
results suggest that the high levels of haplotype sharing 
observed across all nuclear loci are most likely due to 
incomplete lineage sorting and widespread sharing of 
ancestral polymorphism and not due to current gene flow 
between species (Moore 1995; Hoelzer 1997; Omland 
et al. 2006). 

Time since divergence between Icterus 
graduacauda and Icterus chrysater 

According to the IM model, the most probable time of 
divergence between I. graduacauda and I. chrysater was 
estimated to be approximately 300,000 years ago (95% 
HPD 176,000-860,000). This estimated time of divergence 
between these two species is placed during the second 
most recent glacial maximum (Illinoian glaciation), which 
occurred between 130,000 and 300,000 years ago in the 
Pleistocene epoch (1.8 to 0.01 MYA). The Pleistocene was 
characterized by intense climatic oscillations around the 
world, with glacial and interglacial periods that changed 
the distribution of many species (Hewitt 2000). Pollen 
records show that during glacial periods, low tropical for- 
est in southern Mexico were displaced more than once to 
the south during cold periods and to the north during 
warmer periods (Guevara-Chumacero et al. 2010). These 
altitudinal and latitudinal shifts of vegetation were of cru- 
cial importance for avian evolution and speciation, as it 
created the opportunity for species to either move into 
new areas or adapt to changing conditions and in some 
cases go extinct (Hewitt 2000; Zink and Klicka 2006; 
Barrera-Guzman et al. 2011). 
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Dual aspect of the Isthmus of Tehuantepec: 
Barrier and connector 

A variety of studies in phylogeography have focused on 
the Isthmus of Tehuantepec and its influence on the spe- 
ciation of different taxonomic groups including birds 
(Garcfa-Moreno et al. 2004; Cortes-Rodriguez et al. 2008; 
Arbelaez-Cortes et al. 2010; Barrera-Guzman et al. 2011), 
mammals (Sullivan et al. 1997; Kerhoulas 2008; Guevara- 
Chumacero et al. 2010), amphibians (Duellman 1960), 
and plants Palicourea padifolia, (Gutierrez-Rodriguez 
et al. 2011) among others. However, most of these studies 
lacked rigorous multilocus methods to test hypotheses of 
demographic processes (Rosenberg and Nordborg 2002). 
Furthermore, while some of these studies used statistical 
methods such as coalescent analysis, they were carried out 
using mtDNA alone, which can have limitations when 
accounting for stochastic processes and selection, leading 
to different geographic patterns for different loci (Rosen- 
berg and Nordborg 2002; Hickerson et al. 2010; Towes 
and Brelsford 2012). 

Our study shows how utilizing both mtDNA and 
nuDNA offers a more powerful and complete evolution- 
ary history of the species of interest as it compares 
results from both marker types. For example, mtDNA 
sequences showed no intermixing of haplotypes between 
I. graduacauda and I. chrysater, suggesting that the two 
species have been isolated for sufficient time to accumu- 
late fixed mitochondrial differences. Nuclear DNA like- 
wise suggested a long history of divergence with little to 
no gene flow, despite alleles being shared between species 
at all loci. The shared alleles found for nuDNA are likely 
a result of the slower sorting rate of nuDNA relative to 
mtDNA and the retention of ancestral polymorphisms 
(Funk and Omland 2003; Towes and Brelsford 2012). 
The concordant signatures of complete isolation found in 
both marker types strengthen our conclusions about 
history. 

The lowlands of the Isthmus of Tehuantepec have been 
under the influence of climatic changes and shifts of mon- 
tane forests repeatedly over the last million years, creating 
an apparent barrier to gene flow not only for I. graduac- 
auda and I. chrysater but also for other species including 
Chlorospingus ophthalmicus (Bonaccorso et al. 2008), Lam- 
pornis amethystinus (Cortes-Rodriguez et al. 2008), 
Palicourea padifolia (Gutierrez-Rodriguez et al. 2011), 
Peromyscus aztecus (Sullivan et al. 1997) among others. 
However, it has acted as a connector for other species, 
such as: Lepidocolaptes affinis (Arbelaez-Cortes et al. 
2010), Icterus pustulatus, and the Buarremon complex 
(Navarro-Sigiienza et al. 2008). 

In conclusion, the Isthmus of Tehuantepec is the edge 
of the distribution for both species, and there are no 
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historic or present reports of either of them on the other 
side. Not only is the Isthmus an important lowland bar- 
rier, but the narrowing of the land connection forms a 
funnel where only a few individuals close to the range 
edge are near this potential contact zone. The distribu- 
tional allopatry of these two species might be explained by 
different natural history traits, local habitat preferences, 
lack of willingness to cross the lowland habitats of the 
Isthmus of Tehuantepec, competition with resident spe- 
cies, and subtle differences in diet, among others. How- 
ever, it is also possible that interspecific competition for 
resources plays a role. Furthermore, it seems likely that 
differences in plumage coloration may play the crucial 
role in determining whether or not a few individuals that 
cross the lowlands obtain mates and have offspring that 
successfully backcross. Future work could use this interest- 
ing species pair to further explore a range of questions in 
biogeography and speciation. 
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